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Abstract. 

We study the adiabatic quantum dynamics of an anisotropic spin-1 XY chain across 
a second order quantum phase transition. The system is driven out of equilibrium by 
performing a quench on the uniaxial single-spin anisotropy, that is supposed to vary 
linearly in time. We show that, for sufficiently large system sizes, the excess energy after 
the quench admits a non trivial scaling behavior that is not predictable by standard 
Kibble-Zurek arguments for isolated critical points or extended critical regions. This 
emerges from a competing effect of many accessible low-lying excited states, inside the 
whole continuous line of critical points. 
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1. Introduction 

Recent impressive experimental advances in manipulating cold atoms loaded in optical 
lattices have opened up the possibility to investigate the actual dynamics of quantum 
many-body systems with very low dissipation rates and long coherence times PQ; this 
also allowed a very accurate check of the fundamental laws describing the physics of such 
systems. Among the others, it has been possible to probe a variety of very interesting 
and genuinely quantum non-equilibrium phenomena, such as, for example, the collapse 
and the revival of a Bose-Einstein condensate [2], the manipulation of the atomic 
number statistics [3], or the coherent non-equilibrium evolution of one-dimensional 
strongly interacting bosons from a carefully prepared initial state [1]. Furthermore, 
non-equilibrium in cold atomic gases can also be achieved by changing in time some 
of the coupling constants of the system, e.g., the depth of the optical lattice or the 
harmonic trap, on a scale shorter than the relaxation rate. These new experimental 
capabilities have spurred a renewed interest in the study of quantum quenches. 

A lot of attention has been devoted to the study of sudden quenches (see for 
example [5] and references therein). In this paper we deal with an equally debated 
problem, i.e., when the changes in the coupling constants driving the quantum system 
are performed adiabatically. This problem becomes non-trivial if, during the quench, 
the system crosses a Quantum Phase Transition (QPT). Due to the closure of the gap 
in the thermodynamic limit, the system will be unable to stay in its equilibrium ground 
state, no matter how slow is the quench. This problem plays a crucial role in adiabatic 
quantum computation schemes, where the system Hamiltonian is supposed to be slowly 
changed on a time scale that is large as compared to the typical inverse zero-temperature 
gap, so that the system always remains in its instantaneous ground state [El [3, [8]. The 
efficiency of an adiabatic quantum computation algorithm relies on the assumption that 
the minimum gap between the ground state and the first excited state goes gently to zero 
in the thermodynamic limit. When this is not the case, the non-adiabatic evolution close 
to the QPT drives the system out of the ground state. The computation is no-longer 
accurate or, in other words, a number of defects appears in the final state. 

The problem of defect formation in the adiabatic dynamics of critical systems 
was examined much before quantum information: it was first considered by Kibble 
and Zurek (KZ) in the context of phase transitions in the early universe [91 [10] and 
more recently extended to the quantum case [TT1 [12], raising an intense theoretical 
discussion [I3j[ll[ig[lia[171[l8l[ig According 
to the KZ mechanism, the evolution of a quantum system is either adiabatic or impulse, 
depending on the distance from the critical point. The time (i.e., the distance from the 
critical point) at which the system switches from one regime to the other depends on 
the speed of the quench: the slower it is, the later the evolution will become impulse. 
This argument allows to predict the scaling of the density of defects as a function of 
the quench rate. Interestingly, for very slow quenches the quantum evolution can be 
also successfully studied [Tlj [T3] by means of an effective two-level approximation with 



Adiabatic dynamics in a spin-1 chain with uniaxial single-spin anisotropy 



3 



an avoided level crossing, within the Landau-Zener (LZ) formalism [2U E2]. A more 
general scenario arises in the presence of non-isolated quantum critical points, which can 
accumulate and form extended critical regions. Here the validity of the KZ mechanism is 
a priori not obvious, even if in some cases it is still possible to predict the defect density 
by identifying a dominant critical point, or by using scaling arguments [221 E31 EH 125] - 

In this paper we study the adiabatic dynamics in a one-dimensional XY spin-1 
system with single-ion uniaxial anisotropy [331 El] , exhibiting (in equilibrium) a QPT 
of the Berezinskii-Kosterlitz-Thouless (BKT) type. Our interest in the dynamics of 
this specific spin-chain is motivated by the fact that it describes quite accurately the 
properties of the Bose-Hubbard (BH) Hamiltonian both in the limit of strong interaction 
and close to the Mott-to-superfluid QPT [351 [36]; understanding the nonlinear response 
of such system to slow quenches may reveal itself as a powerful tool to probe Bose 
condensates loaded in optical lattices [29]. Interestingly, our results suggest that the 
knowledge of the static properties of the system, in particular of the BKT transition, 
may not be sufficient to predict and fully characterize the dynamical behavior. 

Some dynamical properties of the BH model after a quasi-adiabatic crossing of the 
QPT have been analyzed both from the superfluid to the Mott insulator [30J, and in 
the opposite direction [37], where topological defects arise. Other works focused on 
the emergence of universal dynamical scaling, when quenching to the superfluid phase: 
they started using the original KZ mechanism [121 E]> but then realized that, for non- 
isolated critical points or critical surfaces, a generalization in terms of dynamical critical 
exponents characterizing the whole critical region was necessary [231 Ell ES] . A more 
general analysis of the problem in the context of the breakdown of adiabaticity for 
gapless systems has been presented in Ref. [28]. A numerical analysis of the raising of 
defects in a quenched spin chain model exhibiting a BKT transition has been performed 
in Ref. [22J; in that case defect formation is dominated by an isolated critical point, so 
that a LZ treatment based on the finite-size closure of the dynamical gap at that point 
is still possible. On the other hand, one can also devise a KZ scaling argument, which 
relies on the closing behavior of the gap as a function of the distance from the critical 
point [121 [T7J. I n some circumstances this problem can be quite subtle, since it is possible 
that the gap depends differently on the inverse size of the system and on the parameter 
driving the transition, so that the two approaches give different results: this seems to 
be the case for the system considered in the present paper. We are not aware of further 
quantitative studies of the dynamical defect formation after an adiabatic crossing of the 
BKT transition line; the major obstacle in understanding this type of dynamics raises 
from the fact that here the scaling of defects is generally due to multiple level crossings 
within the whole gapless phase. We believe that this issue deserves further attention. 
This is the aim of the present work. 

The paper is organized as follows. In Sec. [2] we introduce the model and recall the 
main features of its phase diagram. In Sec. [3] we discuss the linear quenching scheme we 
adopt, and define the excess energy of the system with respect to the adiabatic limit: 
this quantity captures the essential physics of the defect formation in the system. All the 
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results of our work are concentrated in Sec. m while in Sec. Owe draw our conclusions. 
2. The Model 

The Bose Hubbard (BH) model [38] , well suited for describing interacting bosons in 
optical lattices [39], is defined by the following Hamiltonian 

^bh = ^2(a\a i+1 +h.c.) + — J^n^rii - 1) . (1) 

i i 

Here a\ (a^) are the boson creation (annihilation) operators on site i (we assumed that 
the lattice is one-dimensional), and rii = a\ai is the corresponding boson occupation 
number. The parameters J and U respectively denote the tunneling between nearest 
neighbor lattice sites and the on-site interaction strength. At integer fillings 1,2,..., 
when the ratio t/U is gradually increased, the BH chain undergoes a QPT of the BKT 
type from a Mott insulating state, where bosons are localized in an incompressible phase, 
to a superfluid, with long range phase order. 

The BH model in equation (TjQ) can be mapped into the effective spin-1 Hamiltonian 
of equation (T2]) in the limit of a large filling and for small particle number fluctuations [35, 
[36]. When number fluctuations are not large it is possible to truncate the local Hilbert 
space to three states with particle numbers n , n ±l (n being the average lattice filling 
per site) . The reduced Hilbert space of site i can then be represented by three commuting 
bosons t a> i {a = —1,0,1), which obey the constraint Yua=-i *a i^«>» = ^ m this wa y' 
the bosons of equation (TjQ) are represented by a\ = y/no + lt\ jio.i + \/^o4,^-i,«- m ^ ne 
limit uq ^> 1 the effective Hamiltonian becomes the Hamiltonian of a one-dimensional 
spin-1 XY chain with single ion anisotropy 



ttspin = -yE W" 5 ^ + h - c -) + D E(^) 2 . ( 2 ) 

i i 

where = \/2(t\ jt 0)i + jt_i,j), = t\ — and with the identification 

^-i, U-D (3) 

(we chose to use the conventional notation for the spin-1 model). In the previous 
equations Sf are spin-1 operators on site i and Sf = Sf ± iSf ; J± and D respectively 
characterize the nearest neighbor coupling strength in the xy plane and an uniaxial 
single-ion anisotropy along the transverse z direction. This system is invariant under 
rotations around the z axis, therefore the total magnetization S* t z ot = £\ Sf is conserved. 
From now on all the quantities are expressed in units of the exchange coupling J± = 1. 

The phase diagram associated to the Hamiltonian in ([2]) is sketched in figure [TJ 
For D > it consists in a large- .D phase for D > D c , that is characterized by zero 
total magnetization (in the limit D —>■ oo each spin has zero magnetization), and a 
BKT transition line for D < D c ; the critical point has been numerically estimated to be 
D c ~ 0.44 [4*01 |4"TI |4"2] . In the rest of the paper we will only concentrate on the adiabatic 
dynamics of the Hamiltonian (T2]). 
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Figure 1. Phase diagram of the spin Hamiltonian @ for J± = 1 and D > 0. 



3. Adiabatic dynamics 

The adiabatic quench is realized by slowly changing the anisotropy parameter D through 
the critical point D c . We suppose to vary D linearly in time: 

D{t) = An " -, With t E [0, r(An - Ain)] ] (4) 
T 

here r is the quenching time scale r, An and -D nn respectively denote the initial and 
the final value of D. In all the cases that will be analyzed we consider An > A, 
and suppose to initialize the system in its ground state; on the other hand we take 
-A n < D c , so that during the quench the system crosses the BKT transition. Since the 
initial ground state has zero total magnetization, and this is conserved by the dynamics 
dictated by equation ([2]), only the excited states carrying zero magnetization will be 
accessible throughout the quench. 

In order to quantify the loss of adiabaticity of the system following the quench, 
we study the behavior of the excess energy with respect to the actual adiabatic ground 
state, after a proper rescaling: 

F (f) (mint) \m) - Waswrnt) h^g) r , 

cxcl ] ~ (H H{t) |Vo> - (1>Qs{t)\H(t) |VM*)> 
where l^o) is the initial state of the system, that is the ground state of Hamiltonian TC(0); 
\ipGs(t)) is the instantaneous ground state of H(t), and \ip(t)) is the instantaneous wave 
function of the system. Strictly speaking, the quantity E eyLC (t) is not defined at the initial 
time t = 0, but one has E exc (t — > + ) = 1; on the other hand at £/ = (An — D^ n )/r, 
the excess energy gives, apart from a constant factor, the final energy cost of defects in 
the system. The final excess energy ranges from E exc (tf) = 1 (totally impulsive case) to 
E e x C {tf) = 0, for a fully adiabatic evolution. 

An exact solution for the spin model in equation ([2]) is not available, not even 
for the static case, therefore one has to resort to numerical techniques. In order to 
investigate both static properties and the dynamics after the quench, we used the 
time-dependent Density Matrix Renormalization Group (t-DMRG) algorithm with open 
boundary conditions [33]. For the dynamics at small sizes L < 10, we checked our t- 
DMRG results with an exact numerical algorithm which does not truncate the Hilbert 
space of the system. For static computations we were able to reach sizes of L = 200, 
while for dynamics simulations we considered systems of up to L = 80 sites. The time 
evolution has been performed with a second order Trotter expansion of H.(t); in most 
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simulations we chose a discretization time step 5t = 10~ 3 , while the truncated Hilbert 
space dimension has been set up to m = 200. 

4. Results 

In this section we describe our results for the adiabatic dynamics of the spin-1 
Hamiltonian. We first analyze the behavior of the excitation gaps which are relevant 
for the quenched dynamics. Then we focus on the dynamics, and discuss the behavior 
of the excess energy as a function of the quenching rate r. We first consider the 
slow-quench region for small system sizes and then concentrate on the scaling regime 
for larger sizes. 

4-1. Dynamical gap 

A great deal of understanding on the adiabatic dynamics derives from the knowledge of 
the finite size scaling of the first excitations gaps. As stated before, since the dynamics 
of the system conserves the total z magnetization, if we suppose to start from the zero- 
magnetization ground state, only excited states with S* ot = will be involved during 
the dynamics. Therefore, the dynamical gap is defined as the first relevant gap for the 
dynamics, that is the energy difference between the ground state and the first excited 
state compatible with the integrals of motion. 

As shown in figure El in the critical region D < D c the dynamical gap A scales 
approximately linearly as a function of the inverse system size L -1 . The same behavior 
also holds for relatively small values of D — D c within the gapped phase, as those 
considered in figure El so that the correlation length is still larger than the size of the 
system, and a quasi-critical regime is found [42J. On the other hand, we numerically 
checked that the leading term in finite-size corrections scale as L~ 2 for D > 2, when 
the system is far from criticality. We extrapolated the value of Ao = A(L — ► oo) in the 
thermodynamic limit by performing a fit of numerical data for L > 50 which includes 
both the leading linear behavior and smaller quadratic corrections. Results are plotted 
in the inset of figure El as a function of D. According to the phase diagram of the 
system, which predicts a closure of the gap for D < D c « 0.44, the asymptotic value 
of the gap is found to be constant and equal to zero for < D < 0.45 (up to values 
~ 2 x 10~ 4 ), while it suddenly raises up as A ~ exp(— c/ \JD — D c ) in the gapped phase 
close to criticality. Thus the dynamical gap closes analogously to the gap between the 
ground state and the first excited state with unconstrained magnetization, which is 
called thermodynamical gap, in a BKT transition [44J. 

The excitation energies of the first three dynamical excited levels in the subspace 
of zero magnetization and for a system of L = 100 sites are displayed in figure [3j as a 
function of the anisotropy D. In the large-/} phase the dynamical gap A^ = A is well 
above the zero; when decreasing D it closes approximately linearly until D ~ 2, then it 
continues closing as far as it approaches a region for D < 0.5, where it becomes almost 
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Figure 2. Ground state excitation energy A in the zero magnetization sector as a 
function of the inverse system size L~ l (we show data ranging from L = 6 to L — 200). 
The various curves are for different values of the single- ion anisotropy D. In the inset 
we plot the asymptotic value Ao in the thermodynamic limit, as extracted from a 
quadratic fit of the data in main panel for L > 50 (black circles) ; the red line displays 
a fit Ao cx exp (-c/y/D - 0.44) of data with c « 2.977. 

constant and very small, as shown in the inset. 

We point out that this type of behavior is quite different from the scenario elucidated 
in the spin-1 /2 Heisenberg model of Ref. [22]. In that case two types of quenches 
involving the antiferromagnetic BKT isotropic point were considered. While in the 
second quenching scheme the system started from the critical region and advanced 
in the opposite direction with respect to our case, the first quench started from 
the antiferromagnetic region and crossed both the BKT point and the ferromagnetic 
isotropic point. Remarkably, the excess energy was found to be essentially characterized 
by the features of the ferromagnetic critical point, where the gap closes faster than in all 
the other points along the critical line. Therefore it was possible to identify a dominant 
critical point which allowed for the applicability of a LZ scaling argument in determining 
the defect density. On the other hand, our quench involves the BKT transition line close 
to the antiferromagnetic isotropic point, and there are no dominant critical points, thus 
leading to a more complex scenario, as explained in the following. 
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Figure 3. Excitation energy A^ of the three lowest excited dynamical levels for 
L = 100 spins in the subspace S* ot = 0, as a function of D; the first excitation energy 
coincides with the dynamical gap: = A. The inset shows a zoom for < D < 1 
of the same plot. 



4-2. Oscillations in the excess energy for slow quenches 

Let us first consider systems of small sizes, as shown in figure H] for L = 6 and L = 8 
sites. We have evaluated the excess energy both with the t-DMRG algorithm (filled 
circles), and with an exact diagonalization which does not truncate the system's Hilbert 
space (empty squares). As the figure shows, data agree well. 

On increasing the rate r, we can recognize two different regimes. For very small 
values of r the excess energy is close to its maximum, and the dependence on the size 
and on r is very small. These points correspond to very fast quenches, where the system 
dynamics is strongly non-adiabatic and the initial state is substantially frozen. As a 
consequence, the state after the quench is found to be in a superposition of many excited 
states of the final Hamiltonian. A second region is characterized by a dominant power- 
law decay, according to E exc ~ r ~ 2 ( see the straight lines in the two insets of figure H]), 
that is superimposed to an oscillatory behavior. This can be explained within a LZ 
approximation: for small values of D and L the gap is large and proportional to 1/L, 
therefore at very small sizes only the ground state and the first excited state participate 
to the evolution of the system, while all the other excited states are not accessible. The 
power-law decay, as well as the oscillations naturally arise when effects of finite duration 
time are taken into account [45]. Following the closing of the gap, the frequency of the 
oscillations decreases at increasing sizes, as it can be seen in the figure. The red curve 
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Figure 4. Final excess energy after an adiabatic quench of D from D m — 1 to 
-Dfin = 0, as a function of the quench velocity r. Filled circles denote t-DMRG data, 
empty squares are obtained with exact diagonalization, while the continuous line is 
a numerical fit with the formula predicted by a LZ model for finite initial and final 
coupling times. The left panel shows data for L = 6 sites, while the right one is for 
L = 8. The insets show the same data in a log- log scale (straight blue lines denote a 
~ r~ 2 behavior). Note the smaller frequency of the oscillations for L = 8. 



displays a fit of numerical data obtained by an effective LZ model in which the initial 
coupling time t{ < is finite, and the final time is tf — (see Appendix A for details 
on the fitting formula). 

The oscillatory behavior can be drastically suppressed starting from a larger value 
of D m , which corresponds, in the LZ model, to decreasing the initial coupling time tf, 
for tf < and in the limit of £, — > — oo the oscillations disappear and a pure power-law 
~ t -2 decay survives [45]. This is seen to emerge from numerical data of figure EJ where 
we started quenching from D m = 4. Notice also the substantial independence of E exc 
on the system size in the fast quenching limit. 

4-3. Scaling regime 

The analysis of the effects of the quantum phase transitions on the adiabatic quench 
dynamics demands sufficiently large system sizes. We now concentrate on this aspect 
and study the excess energy as a function of r for considerably larger values of L. Due to 
the increasing computational difficulty in simulating large systems, we restrict ourselves 
to quenching schemes in which = 1. 
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Figure 5. Excess energy after an adiabatic quench of D from D m = 4 to _Dg n = 0, as 
a function of r. The various data are for different system sizes L — 6 (black circles), 
8 (red squares), 12 (green diamonds). The straight blue line indicates a behavior 
E eX c ~ t~ 2 and is plotted as a guideline. 



In figure [6] we plot the final excess energy of the system after a quench from D m = 1 
to Df[ n = of time duration r. Starting from fast quenches and going towards slower 
ones, we can now distinguish three different regimes: the first strongly non-adiabatic 
regime at small r is analogous to the one previously discussed for small sizes. In the 
opposite limit of very slow quenches r ^> 1, we also recover the power-law r~ 2 behavior 
superimposed to oscillations coming from an effective LZ description with finite coupling 
duration. Most interestingly, in between these two opposite situations, a characteristic 
power-law regime emerges, where: 

£cxc~r- Q with ae [1,2]. (6) 

This is dominated by transitions to the lowest dynamically accessible gap, and it is 
crucially affected by the critical properties of the system. The crossover time r* at 
which this regime ends typically increases with the size, as it can be qualitatively seen 
from the figure (arrows denote a rough estimate of r* for the different sizes), and diverges 
in the thermodynamic limit; unfortunately we were not able to analyze the scaling with 
L, because of the intrinsic difficulty in estimating the ending point of the r~ a behavior. 
Nonetheless, even at asymptotically small quenching velocities, for very large sizes the 
scaling of defects (El) ruled by criticality persists, thus meaning that the system dynamics 
cannot be strictly adiabatic. 

The scaling of the decay rate a with the size has been analyzed numerically, for 




Figure 6. Final excess energy as a function of the quench rate r. The various panels 
stand for different system sizes. Symbols denote numerical t-DMRG data, while the 
straight line is a power-law fit that has been performed for r < r* (r* is indicated by 
the vertical arrow). The two straight dashed lines in the upper panels denote a r~ 2 
behavior, and are plotted as guidelines. The oscillating dashed line at L — 10 is a fit of 
data with big r, according to the LZ model for a finite coupling duration. The values 
of a corresponding to the power-law fits are quoted in each panel. Here we set D m = 1 
and -Dfin = 0. 



data corresponding to L ranging from 10 to 60 spins; at L < 10 this regime was not 
identifiable. Some representative cases are shown in figure where each of the four 
panels stands for a given system size, while straight continuous lines indicate the best 
power-law fits of the scaling regions. In the case of L = 10 sites (upper left panel), we 
cannot give a reliable estimate of a, since the width of the scaling region is narrow and 
the fit is very sensitive to its actual starting and ending points. The straight line in the 
plot corresponds to a ~ 1.798 and has been obtained from a power-law fit of numerical 
data from r = 1 to r* = 3. As one can see, this is hardly distinguishable from the t~ 2 
power-law behavior of the slow-quench regime (straight dashed line), thus meaning that 
the existence of the scaling region itself is here in doubt. This is not the case for the 
other panels, where a power-law behavior of the type in equation is clearly visible. 
Namely, we fitted our data until the r* value, that is labeled in figure [6] by a vertical 
arrow: as we could expect, the size of the scaling region increases with L. 

Summarizing the results obtained for the various sizes, in figure [7] we report the 
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Figure 7. Power-law decay rate a in the intermediate scaling region for the excess 
energy, as a function of the system size L. The quench is performed from D m = 1 to 
£*fin = 0. In the inset we plot the same quantity as a function of 1/L. The blue line is 
a linear fit of data with L > 30, and predicts an asymptotic value of s» 1.28 in the 
thermodynamic limit. 



behavior of a as a function of L (in the inset we plot the same data with 1/L on the 
x-axis). The uncertainty affecting the value of a extracted from the power- law fits of 
numerical t-DMRG data is mostly due to the inaccurate knowledge of the extremes of 
the scaling region. For each value of L, we identified a trial power-law region and then 
computed several values of a by progressively sweeping out the points from that region, 
starting from the borders. We then evaluated error bars, that are displayed in the plot, 
by performing a statistical analysis of the values of a thus obtained. In order to give an 
estimate of the power-law decay rate in the thermodynamic limit, we supposed that, at 
large L, a scales inversely proportional with the system size. In this way, performing 
a linear fit of data with L > 30, we extracted the asymptotic value ~ 1.28 in the 
thermodynamic limit (see straight blue line in the inset). 

We would like to stress that, in this context, our numerical results seem not to find 
a straightforward explanation with LZ or KZ arguments. One could, for example, try 
to follow a standard LZ argument, that is based on the assumptions that the dynamical 
gap scales linearly with the inverse size, nearby and inside the critical region, and that 
the adiabaticity loss is essentially due to the presence of a dominant critical point, where 
the gap closes faster than elsewhere [11], [22] . In the LZ approximation, the probability 
of exciting the ground state is a global function of the product tA^, where A m is the 
minimum gap achieved by the system during the quench. Assuming a critical scaling 
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Figure 8. Power-law decay rate a as a function of the ending value for the quench 
.Dfin < D c and for fixed D m = 1. Data are for different system sizes, as explained in 
the caption. In the inset we display the excess energy for a quench ending at various 
Dfi n , and for a system size L = 60. straight lines are power-law fits in the scaling 
regime. 



of the gap A m ~ as shown in figure El the density of defects can be estimated by 
evaluating the typical length L E of a defect-free region, once the probability for this to 
occur is e. As a consequence, this would give E exc ~ l/L £ oc r -1 / 2 , exactly as in the Ising 
case, in contrast with numerical evidence. On the other hand, a scaling argument based 
on the KZ mechanism can be adopted [TTJ [12] ; this relies on the fact that the dynamical 
gap A seems to close with D — Dc as the thermodynamical gap in a BKT transition, 
that is, it depends on the anisotropy parameter as A oc exp(— c/y/D — Dc) |44j, 
so that the critical exponent for the correlation length v diverges. In this case, the 
KZ scaling argument predicts a power-law scaling exponent a = (d + z)vj{zv + 1), 
d being the dimension of the system and z, v critical exponents; in our specific case 
d = z = l,v — > oo thus leading to a = 2 (plus some logarithmic corrections) [12j [T7] . 
This again contrasts with our numerical evidence, thus revealing that the presence of 
a critical line in which the gap closes always in the same way seems to indicate that 
all the low-lying excitation spectrum becomes necessary to predict the actual behavior. 
We notice that the two above mentioned different dependencies of the dynamical gap 
on the inverse size, like ~ 1/L for L — > oo, and on the distance from the critical point, 
as ~ e ~ c /VD-D c ^ are con fi rmec i q U ite precisely by our data. 

We checked the dependence of a on the final point of the quench: in figure [8] we 
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varied the ending point D &n , while keeping D in and the system size fixed (explicit data 
for the excess energy E exc as a function of r are presented in the inset, at L = 60). For 
values of D^ n outside the critical region we find that a depends on D^ n at finite sizes. 
Nevertheless, the range of the scaling region shrinks with L and eventually disappears 
in the thermodynamic limit, so we argue that the dependence of a on D^ n > D c should 
be entirely due to finite size effects. For D &n < D c we observe that the dependence 
of a on Dg n weakens as the system size is increased. In this case the scaling region is 
valid until a quench rate r* +00; the power-law decay rate tends to a value that is 
independent of D^ n < D c and has been extrapolated from numerical data of figure [7] to 
be doc ps 1.28. 

5. Conclusions 

In this work we have analyzed the quenched dynamics of a quantum anisotropic spin-1 
XY chain, when it crosses a Berezinskii-Kosterlitz-Thouless quantum phase transition. 
The quench has been performed on the uniaxial single-spin anisotropy, and has been 
chosen to vary linearly in time with a given velocity. We focused on the residual excess 
energy of the system after the quench, and studied its dependence on the velocity of 
the quench. For very slow quenches and finite system sizes we were able to describe 
the properties of the system in terms of an effective Landau Zener model, where the 
system can only get excited to its first excited state. Most interestingly, we pointed out 
the emergence of an intermediate region where the excess energy drops as a power-law 
with the quench rate, and exhibits a non trivial scaling behavior. At least for the finite 
sizes considered here, the decay rate depends on the size of the crossed critical region, 
and cannot be explained in terms of usual scaling arguments, such as the standard 
Kibble-Zurek mechanism and its generalization to critical surfaces [221 1231 [241 125] . In 
the thermodynamic limit the system obeys a non-trivial scaling behavior E exc ~ t~ q , 
with 1 < a < 2, even when r — > 00 (i.e., for very slow quenches). 
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Appendix A. Landau-Zener model for finite coupling times 

The Landau-Zener (LZ) model consists in a two level system describing an avoided 
level crossing: two energy levels moving in time are widely separated at first, then they 
approach each other with time, and finally part away again [3UI32]. When the two levels 
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are well separated, each eigenstate preserves an individual character; on the other hand, 
when levels are close together, they mix due to their interaction. The Hamiltonian is 
given by: 

with a detuning A(t) = f3 2 t (where [3 2 > 0), and a time independent coupling fl that, 
in the original LZ model is supposed to last from ti = — oo to tf = +00 [311 [32] . Here 
we review the general case where the coupling is turned on at ti and off at tf [45J. 
The equation flA.lj) is written in the basis of the two eigenstates of the Hamiltonian in 
absence of interaction. 

The probability amplitudes C{tA = [Ci(ij), C 2 (ti)] T for the two levels at the 
beginning are connected to the ones at the final time tf by the unitary evolution matrix 
U(tf,ti), so that: C(tf) = U(tf,ti) C{tj). Their elements are given by: 

lYi - Urn 2 ) 

U 11 (T f ,T i ) = V 2 V ^M^e--/ 4 ) (A.2) 

x D_ 1+luJ 2 /2 (T t V2e-™^)} , 
U 12 (T f ,Ti) = F(1 - Jl^ e^ 4 [-D^ /2 (TfV2e-^) (A.3) 



+ D^ /2 (T f V2e t3n/4 ) D^/^T^e-^)} , 

where we have introduced the rescaled time T = j3t and the scaled dimensionless 
coupling strength uj = while D v (z) denote the parabolic cylinder functions. If 

we suppose that the system is initialized in its ground state, i.e., Ci(tj) = 1, C^ifi) = 0, 
the transition probability to the excited state at the final time is given by pW(t f ,ti) = 
\U 2 i{tf,U)\ 2 . 

This is related to the relevant adiabatic basis, that is the basis of the instantaneous 
system eigenstates, by a unitary transformation. If A(t) = [Ai(t), A 2 (t)] T are the 
probability amplitudes for the two levels in the adiabatic basis, then A(t) = R(t) C(t), 
where R(t) is the rotation matrix 

R(T) -|^sin^) cos^(t) J' (A ' 4) 

with tan[2i?(t)] = Q(t)/A(t). Therefore, the evolution matrix in the adiabatic 
representation is given by U a (tf,ti) = H T (tf) U(tf,ti) R(t;), and the adiabatic-following 
solution for the transition probability is P( a \tf,ti) = \U^ (tf,ti)\ 2 . 

For the original LZ model, where the coupling is supposed to last from ti — > — 00 
to tf — > +00, the expression for the excitation probability at the end of the quench in 
the adiabatic basis simplifies to an exponential form: 

P (a) (+oo,-oo) =e-™\ (A.5) 
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In the case of a finite coupling duration, that ends before or exactly at the crossing 
(i.e., tf < 0), we have a much involved expression, which predicts a leading power-law 
behavior ~ r~ 2 superimposed to an oscillating behavior. Eventually oscillations 
are damped for long lasting couplings: in the limiting case where the quench ends at 
the critical point and is infinite lasting (tj = — oo, tf = 0), the probability is given by 

pla ^-^ = Tb'^- (A - 6) 

The scaling with the quench velocity r follows from the fact that the times t oc r, while 
f3 2 oc 1/t (this implies that u oc \fr). 

We used the explicit formula for the adiabatic transition probability P^ a \0,ti < 0) 
in order to fit t-DMRG data for the excess energy of our system in the regime of large 
r, where defects still do not form and the quench dynamics can be considered adiabatic. 
While it is clear that i; < in our case, it is not obvious a priori whether tf < or 
tf = 0, since LZ relies on the assumption that there is only one point of closest approach 
of the energy levels; on the contrary, in our model we have a whole critical line. We 
actually chose tf = and used tj < as a fitting parameter, having not a rigorous 
criterion at our disposal, but following the qualitative picture suggested from figure [3j 
the gap closes monotonically during the quench, reaching the minimum at D^ n . The 
red curves in figure H] have been obtained by fitting numerical data with the theoretical 
prediction given by P'"'; we admitted a global rescaling pref actor qb and imposed the 
following constraints: T, = —T y/r, Tf — 0, uj — u y/r. The fitting parameters are 
T , luq, (f). For the left panel (N = 6) we chose T « 0.8, Uq = 0.84, qb = 1.72, while for 
the right one (N = 8) T « 0.735, u = 0.7, = 0.9. 
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